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PREFACE 


Part  of  the  Project  RAND  research  program  consists  of 
basic  supporting  studies  in  mathematics.  This  Memorandum 
presents  a  method  for  the  computational  solution  of  a 
type  of  differential  equation  that  frequently  arises  in 
the  construction  of  mathematical  models  in  a  variety  of 
fields,  including  engineering  and  biophysics. 


SUMMARY 


Functional  differential  equations  of  the  form 

(1)  u'(t)  -  g(t,u(t),u(h(t))), 

and,  more  generally,  of  the  form 

(2)  u'(t)  -  g(t,u(t),u(h(u,t))), 

arise  in  the  construction  of  realistic  models  in  a  number 
of  fields,  ranging  from  electromagnetic  theory  and  control 
theory  to  respiratory  theory  and  neurophysiology.  The 
analytic  aspects  are  quite  complex,  and  numerical  solution 
Is  anything  but  routine,  even  with  the  aid  of  a  digital 
computer.  In  this  paper,  we  wish  to  extend  a  method  for 
the  conp  utational  treatment  of  differential~dif ference 
equations  to  cover  equations  of  the  form  given  in  (1). 
Equations  of  the  type  appearing  in  (2)  can  then  be 
treated  by  means  of  successive  approximations. 
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ON  THE  COMPUTATIONAL  SOLUTION  OF  A  CLASS  OF 
FUNCTIONAL  DIFFERENTIAL  EQUATIONS 

1 .  INTRODUCTION 

Functional  differential  equations  of  the  form 

(1.1)  u'(t)  -  g(t,u(t),u(h(t))), 

and,  more  generally,  of  the  form 

(1.2)  u' (t)  -  g(t,u(t),u(h(u,t))), 

arise  in  the  construction  of  realistic  models  in  a  number 
of  fields,  ranging  from  electromagnetic  theory  and  control 
theory  to  respiratory  theory  and  neurophysiology.  The 
analytic  aspects  are  quite  complex  (see  [1],  [2]),  and 
numerical  solution  is  anything  but  routine,  even  with  the 
aid  of  a  digital  computer.  In  this  paper,  we  wish  to 
extend  a  method  given  in  [3,4,5]  for  the  computational 
treatment  of  differential— difference  equations,  to  cover 
equations  of  the  form  given  in  (1.1).  Equations  of  the 
type  appearing  in  (1.2)  can  then  be  treated  by  means  of 
successive  approximations. 

2^l  EQUATION  u'(t)  »  g(u(t)  ,u(h(t))) 

Let  us  suppose  that  h(t)  ^  t,  for  t  that 

the  future  is  determined  by  the  past.  An  objection  to 
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uslng  straightforward  techniques  for  the  numerical 
solution  of 

(2.1)  u'(t)  -  g(u(t),u(h(t))), 

with  appropriate  initial  conditions,  lies  in  the  amount 
of  rapid-access  storage  required  for  large  t.  If  this 
equation  is  merely  one  of  a  number  of  equations  which 
arise  in  the  description  of  a  large,  complex  system,  we 
cannot  afford  to  use  up  rapid-access  storage. 

Let  us  suppose  that  h'(t)  •  0  for  t  _  0,  and 

let  the  inverse  function  h"~^(t)  be  denoted  by  H(t), 
Suppose  further  that  u(t)  is  known  in  some  initial 
interval  [0,t^],  where  tj^  «  H(0),  and  let  the  sequence 
[t^]  be  defined  recursively  by 

(2.2)  t^  -  H(t^_j^),  n  -  2,3,.... 

Let  H^*^^(t)  denote  the  n— th  iterate  of  H(t), 

n  =  2,3,....  We  observe  that  the 
function  H(t)  maps  the  inteip^al  onto 

^^n^^n+L^  in  one— to— one  fashion,  n  =  1,2,...,  and 

maps  onto  I  ’ 

Consider  the  function 

=  u(M^'’^(s)),  n  =  0,1,2,..., 


(2.3) 


where  s  is  restricted  to  the  interval  [0,tj^],  and 
H^^^(s)  ■  s.  Thus  the  values  of  u^(s),  for  0  <  s  <  tj^, 
are  the  values  of  u(t)  for  t^  1  1  ‘n+r  We  have 

(2.4)  u^Cs)  -  H<'’)(s)]u'(H<">(s)), 

and  the  derivative  of  H'^^(s)  can  easily  be  evaluated 
recursively  by  the  formula 

(2.5)  h(")(s)  -  H'(h("-^)(s))  4  H<"“^^(s). 

Now  we  set  t  •  H^'^^(s),  where  0  <  s  <  t..  Then 
from  the  equation  in  (2.1)  we  get 

(2.6)  u'(H^^^(s))  -  g{u(H^^^  (s)),u(H^^^^  (s))} 

■  g(Uj^(s),u^_3^(s);,  n  -  1,2,.... 
Referring  to  (2.4),  we  see  that  (2.6)  may  be  written 

(2.7)  (s)]g(u^(s),u^j^(s)),  n  -  1,2,... 

Thus  (2.1)  has  been  replaced  by  a  system  of  ordinary 
differential  equations  where  s  now  ranges  over  a  fixed 
interval  [0,t^].  This  is  important  from  the  computational 
point  of  view. 


3.  NUMERICAL  SOLUTION 


The  computational  solution  of  this  system  is  not 
routine,  since  we  do  not  possess  the  requisite  initial 
values  [u^(0)}  for  n  ^  2.  Consequently,  we  proceed  as 
in  [3].  The  equation 

(3.1)  u|(s)  =  H' (s)g(u^(s),UQ(s)),  u^(0)  =  UQ(t^) 

gives  us  the  value  We  now  consider  the 

simultaneous  equations 

(3.2)  uj(s)  «=  H' (s)g(uj^(s),Uy(s)),  Uj^(O)  =  UQ(tj^), 

u^(s)  =  (s)] 'g(u2(s),u^(s)),  U2(0)  =u^(t 

and  solve  these  to  obtain  u^(0)  =  U2(tj^), 

Continuing  in  this  way,  we  require  no  storage  of 
functions,  at  the  expense  of  being  required  to  solve 
successively  larger  systems  of  equations, 

4.  ANOTHER J^TOOD 

Reduction  of  equation  (2.1)  to  a  system  of  ordinary 
differential  equations  over  a  fixed  interval  can  also  be 
achieved  by  introducing  new  independent  and  dependent 
variables  in  a  manner  attributed  originally  to  Laplace 
(see  [1],  page  84).  It  is  also  closely  related  to  some 


ideas  of  Abel  concerning  the  iteration  of  functions.  We 
begin  by  introducing  a  new  independent  variable  x  by 
the  equations 

(4.1)  t  =  p(x),  h(t)  =  p(x  -  §), 

where  ?  is  an  arbitrary  fixed  constant.  The  function 
p  must  be  chosen  so  that 

(4.2)  h(p(x))  «=  p(x  -  ?), 

or 

(4.3)  p(x)  =  H(p(x  -  §)), 

and  should  be  monotone  and  continuous,  in  order  that  the 
relation  t  =  p(x)  be  solvable  for  x  =  P~^(t),  We  note 
that  if  p  is  any  monotone  function  on  0  <  x  <  '  with 
the  property  that 

(4.4)  h(p(‘"))  =  p(0), 

then  because  H  is  continuous  and  monotone,  it  will 
follow  that  p(x)  is  continuous  and  monotone  on  0  <  x  <  2 
and  by  iteration  of  (4.3),  continuous  and  monotone  on 


Let 


(4.5)  v(x)  -  u(t)  -  u(p(x)). 

It  then  follows  that 

v'(x)  -  u'(t)p'(x)  -  u' (p(x))p'(x). 

Since  (2.1)  yields 

u'(p(x))  -  g(u(p(x)),u(h(p(x)))) 

-  g[u(p(x)),u(p(x  -  0)}# 

and  since  u(p(x))  ■  v(x)  and  u(p(x  ~  ?))  ■  v(x  —  ^), 
we  see  that 

(4.6)  v'(x)  »  p' (x)g(v(x),v(x  -  ?)). 

Thus  the  introduction  of  new  independent  and  dependent 
variables  x  and  v  by  means  of  (4.1)  and  (4.5)  leads 
to  the  replacement  of  (2.1)  by  the  pair  of  equations 
(4.2)  and  (4.6),  one  of  which  is  a  difference  equation 
and  the  other  a  differential— difference  equation,  both 
with  a  fixed  lag  %  It  is  not  hard  to  see  that  if  p 
and  V  are  continuous  solutions  of  (4.2)  and  (4.6),  then 
(4.1)  and  (4.5)  define  a  solution  u(t)  of  (2.1). 
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The  equations  (4.2)  and  (4.6)  can  now  be  reduced  to 
a  system  of  equations  on  the  interval  (0,?]  by  the 
technique  of  [4].  That  is,  we  define 

(4.7)  “  P^^ 

(0  ^  X  <  n  *  0,1,...) 

V  (x)  “v(x  +  nO. 
n 

Here  by  Pq(x)  and  Vq(x)  we  accordingly  mean  the 
initial  values  of  p(x)  and  v(x)  on  0  <  x  <  ?.  With 
these  definitions,  equations  (4.3)  and  (4.6)  finally  take 
the  form 

(4.8)  Pn^^^  '  <  ^  <  “  1,2,...) 

(4.9)  “  Pn^^)8(Vj^(x),v^_j^(x)). 

The  solution  of  (4.8)  can  now  be  carried  out  by  iteration, 
and  the  solution  of  (4.9)  by  the  method  sketched  in 
Sec.  3.  The  function  P^C^)  evaluated  recursively 

from  the  formula 

p^(x)  -  H' (p^_j^(x))p^_j^(x). 

The  similarity  of  (4.8)  to  the  relation 
H^^^(s)  =  (s)) ,  and  of  (4.9)  to  (2.7),  is 


* 

apparent.  Indeed,  If  we  choose  Pq  so  that  Pq(x)  ■  x, 
then  from  (4.4)  we  obtain  h(f)  ■  0,  and  ^  Is  the  same 
as  the  tj^  defined  in  Sec.  2.  Moreover,  we  then  obtain 
p^(x)  ■  and  since 

v^(x)  -  v(x  +  n?)  -  u(p(x  +  nO)  "  u(p^(x)), 

the  function  v  (x)  is  identical  to  the  function  u  (s) 

n  n ' 

of  Sec.  2. 

On  the  other  hand,  the  formulation  in  this  section 
seemingly  allows  some  extra  latitude,  since  the  choice  of 
Pq(x)  is  largely  arbitrary.  In  practice,  however,  it  is 
unlikely  that  any  advantage  can  be  derived  from  this, 

5.  AN  EXAMPLE 

We  wish  to  consider  a  particular  case,  in  order  to 
illustrate  a  possible  pitfall  in  the  application  of  the 
techniques  given  above.  Suppose  the  equation  in  question 
is 

(5.1)  u'(t)  -  g(u(t),u(-|  -  1)), 

so  that  h(t)  *  t/2  —  1,  H(s)  «  2(s  +  1).  Then  t^^  * 
t2  “6,  «  14,  and  so  on,  and 

H^"^(s)  -  2"s  +  t^. 


2, 
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Equation  (2.7)  becomes 

Uj^(s)  -  2"g(u^(s),u^j^(s)). 

Now  suppose  we  carry  out  the  solution  of  this  system, 

determining  equally  spaced  points 

s  «  where  N6  -  tj^.  Referring  to  (2.3),  we 

see  that  we  obtain  in  this  way  the  values  of  u(t)  at 

the  points  t  ,t  +  2'^fi,t  +  2*^2fi,..,,t  ...  Since  the 
^  n  n  n  n+1 

length  of  the  Interval  ^^n^^n+1^  doubles  when  n 
increases  by  one,  whereas  the  number  of  points  at  which 
we  know  u(t)  is  unchanged,  it  is  evident  that  our 
computed  values  provide  less  and  less  information  about 
u(t)  the  larger  t  becomes. 

Two  ways  of  overcoming  this  difficulty  suggest 
themselves.  One  is  simply  to  use  a  very  small  value  of 
4  in  the  first  place — but  this  exacts  a  penalty  in 
storage  and  computing  time.  The  second  is  to  combine  the 
above  scheme  with  an  interpolation  process.  For  the 
equation  in  (5.1),  this  could  take  the  following  form. 
Starting  with  Uq(s),  tabulated  at  0,fi,...,N6,  solve 
the  differential  equation  to  obtain  Uj^(s),  and 
consequently  u(t)  at  N  points  on  the  interval  [tj^,t2l. 
By  interpolation,  subtabulate  u(t)  at  N  additional 
points.  This  yields  Uj^(s)  at  2N  points,  say 
“  0, fi/2, 4, 36/2, . . . ,N6.  Continuing  in  this  way,  we 
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obtain  u(t)  at  equally  spaced  points  over  any  interval. 
It  seems  likely  that  the  insertion  of  an  interpolation  at 
each  step  of  the  process,  or  whenever  the  computed  points 
become  too  thinly  distributed,  will  improve  the  accuracy 
of  the  results. 

In  many  applications  of  interest,  fortunately, 
h(t)  is  such  that  h(t)  ~  t  —  b,  b  a  constant,  as 
t  -•  OD  .  Hence,  the  foregoing  difficulty  does  not  arise. 

6 .  SUCGE SSIVE  APPROX IK^^^ 

A  minor  modification  enables  us  to  treat  the  more 
general  equation  of  (1.1).  The  equation  of  (1.2),  which 
arises  from  some  realistic  models  of  respiratory 
processes,  can  be  treated  by  means  of  successive 
approximations . 
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